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INRODUCTION 

Nowadays Molecular Dynamics (MD) approach is 
widely applied in many areas of physics, chemistry and 
biochemistry. The MD is based on solution of second or- 
der differential equations of motion, this is why the integra- 
tion algorithm is a cornerstone of the MD method 1 1 1. The 
Newtonian equations of motion are time-reversible and it 
would be reasonable to preserve this essential property in 
our integration schemes. Since 1990 there are many nice 
symplectic integrators were invented |2|, mainly in force- 
gradient form but none with higher-order gradient opera- 
tors. In this short report we derive a new time-reversible 
explicit integrator on the basis of second order Tailor ex- 
pansion of force. There is good reason to think the new 
method will be easy-to-use for MD and, possibly, celestial 
mechanics applications. 



IDEA OF A NEW METHOD 

Consider the second order differential equation: 

x = f(x) 

It is useful to introduce the function s(t) : 

1 + t/h, -h < t < 



(1) 



8 ® = \l-t/h, 0<t<h & 
By using Eq.Q one may integrate by parts the integral 



h p0 fh 

f(x)s(t)dt= / f(x)s{t)dt+ / f(x)s(t)dt 

-h J-h JO 

= (xs — xs)\°_ h + (xs — xs)\q, 



and finally 

x(-h) - 2x(0) + x(h) = h [ f(x)s(t)dt. (3) 

J-h 

A proper approximation of the function f(x) in the Eq.Q 
within the segment t € [—h, h] may give us difference 
schemes for numerical integration of Eq.Q. For instance, 
by assuming f(x) « f(x(0)) , the exact formula Q im- 
mediately gives the explicit Verlet integrator: 



x(-h) - 2x(0) + x(h) = f(x(0))h 2 



(4) 



Because Eq.@ is symmetric for t ± h the Verlet method 
is time-reversible. In the case of time quadratic approxima- 
tion of f(x(t)) the implicit Stoermer method is derived 
from Eq.Q: 



x(-h) - 2x(0)+x(h) = 
= (f(x(-hj) + 10/(a;(0)) + f(x(h)))h 2 / 12 



It is convenient to introduce notations: 

x = x(Q), x h = x(h), x-h = x(—h) 

5x = x - x-h, Sx h =Xh—Xo 

v(t) = x(t), a(t) = x(t), b(t) = a{t) (6) 

fo = f(xo), fh = f(x h ), f-h = f(x-h) 

f'(x) - df(x)/dx, f"{x) = d 2 f(x)/dx 2 



Thus the Verlet method @ is given by 

Sxh — 5xq + aah 2 + 0(/i 4 ) 



(7) 



Expand f(x) in the vicinity of x — xq and x(t) 
at a point t = in a Taylor series 

f(x) = fo+fo(x-x )+f '(x-xo) 2 /2+f "(x-xof/6+.. 

(8) 



x - x Q = v a t + a Q t 2 /2 + b Q t 3 /6 + .. 



(9) 



Substitute (x — xq) from Eq.(|9} to Eq.(|8j and hold 
only the even terms so one may find that the function f even 
along the trajectory of motion around t = is 

feven(t) = f + (<*o/o + iftfg)* /* + O(^) (10) 

Substitute Ea.dlOl to Eq.0 and integrate the latter: 

r h 

h I t 2 s(t)dt = h 4 /6 

-h 



and 



co = a /o + Vofd 



S Xh = 6x + a h 2 + c Q h 4 /12 + 0(h 6 ) (11) 

The explicit integrator il It is time-reversible likewise the 
Verlet method. The important difference of Ea. il 1> from 
Verlet formula is a velocity and acceleration dependen- 
cies in the cq coefficient. 

How to evaluate velocity Vh ? 

1st way. By using the same approach as above for 
derivation of ( fTTl one may deduce a formula 



Vh = v-h + I f(x)dt. 
'-h 



(12) 



(5) 



Therefore the time-reversible velocity formula is 

v h = v- h + f 2h + (a /o + v 2 fl!)h 3 /3 + 0{h 5 ) (13) 

The main disadvantage of Ea. (ll3> is a poor accuracy. 

2nd way. Assume that in the vicinity of t = the 
x(t) can be represented by polynomial interpolation: 

x(t) =x Q + v t + a Q t 2 /2 + c 3 t 3 + c 4 t 4 + c 5 t 5 + c 6 t 6 

By using the known positions, velocities, and accelerations 
at points t = [— h,0, h] (see Table 1) one can derive a 



Table 1 The source data for the time-reversible interpola- 
tion formula of velocity (II 41 
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Table 2 Algorithm of integrator based on formulae (II 1I15> . 
Here K is a time-step number and L is a logical step 
number within the given time-step. The first time-step 
Oq — > 1q has to be done by using another integrator. 
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time-reversible formula for evaluation of velocity at time 
point t = h : 

8 h 
Vh = V- h + (5x h -5x )-^ + (a- h + a h -8a Q )- + O(h (i ) 

(14) 

Ea. (ll4> is more accurate than Eq. (ll3> . but it requires the 
knowledge of position and acceleration at t = h . The 
general scheme of the new integrator on the basis of formu- 
lae (II 11141 is presented in the Table 2. The first time-step 
has to be done by using another integrator. 

3rd way. One may obtain even more precise velocity 
formula by using correction term cq of Eg. (II 1>: 



Vh = V-h + (Sx h - fix )j— + (a- h + a h 



24a )A- 



Co 



8h 3 
12 • 13 



0(h 8 ) 



(15) 



It should be noted the precision of estimate of velocity 
by Eq.(ll5> exceeds the precision of coordinate evaluation, 
therefore the energy conservation can not be improved con- 
siderably by this way. 

VECTOR FORMULAE 

Again let us consider the second order differential 
equation in vector notation: 



r = f (r) /m = a 
Eg s . 1218191 can be rewritten as vector equations: 



(16) 



As it was for Eq. fllOt . substitute Ar from Eq.(ll9> to 
Ea.(ll8> and hold only the even terms 



f(r + Ar) even = f + [(ao • V) f + (v • V) 2 f ] 

+ 0{t A ) 



and finally 



c = [(a-V)f+(v.V) 2 f] /m 
5r h = Sr Q + aah 2 + cq^/U + 0{h 6 



2 + 
(20) 

(21) 



For evaluation of velocity one may again use Eq.( ll5l l and 
the general scheme of integration (see Table 2). 
1 st example 

Consider the forces depending only on distance r = |r| 

f(r) - f(r) r 
We need the following vector identities 

f 

(w-V)f = f(r) w+— (w-r)r (22) 
r 



f f 
(w • V) 2 f = 2 — (w • r) w + — (w • w) r- 

r r 

f" f \ , N 2 



(23) 



where w = const . 

Applying Eqs.(|22j23} to the general integration for- 
mula Eq. (l21> 



c 



/(r)a + 2-(vr)v+^(v 2 - 

r r 



+ (a.r))r+(4-4)(v-r) 2 r 



(24) 



/m 



5r h = Sr + a h 2 + c /i 4 /12 (25) 
2nd example 

Let us consider the Molecular Dynamics system of N par- 
ticles with pair-wise interaction. In the vector notation N 
equations of motion are given by 



(26) 



Represent the relative positions in a Taylor series 
about t = 



h f h 

Sr h = Sr + — / f(r)s(t)dt. (17) Ar^ = ry(t)-r y (0) = (vi-Vj) i+(a l -a i ) t 2 /2 + 

m J-h (T 



f(r + Ar) = f + (Ar-V) f + (Ar-V) 2 f /2 + ... (18) 



Ar = v t + a < 2 /2 + b a t 3 /6 + 



(19) 



(27) 

and expand the pair-wise forces fy in the vicinity of 

r«(0) 

(0) + A 

+ (Ar • V) 2 ? % /2 + .. 
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exponent n in timestep h = 2 ~ n , -\0Q 2 h 



Fig. 1 Time-averaging numerical accuracy of two methods 
as a function of time-step at equal numbers of MD simula- 
tion steps. 
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Fig. 2 Comparison of the relative CPU times r required 
for simulation of 1 MD time unit at given accuracy. It is 
assumed the r = 1 at time-step h = 1 for Verlet method 
and the r = 2 at time-step h = 1 for a new method. 



By substitution Ary from Ea. J27t to Ea.J28> one may 
obtain 



Sr ih = 8r i0 + V* — f, 

f * m - 



h 4 



ij,0 



12 m,- 



:ifi • - ' (29) 

((a» - a,-) ■ Vy) fy,o + (( v ; - v j) • Vy)o fy,o] 

Then, by substitution ay = (a, — a,) and so on into 
Eg s. ( 1241231 it is easy to derive the final formula: 



Cy,rj - 



f 

f(r) a +2— (vr) v- 
r 



+ ^(v 2 + (a-r))r- 

r 



rf'-f' , 2 



-(v • r) r 



I mi 



ij,0 



(30) 



5r irh = 5r hQ + a hQ h 2 + ^3,0 h 4 /l2 (31) 

MOLECULAR DYNAMICS TESTING 

We have performed a practical test of the new method 
to estimate its applicability for MD simulation problem. 
The pair potential for our model system is described by the 
quasi-Lennard-Jones (qLJ) pair potential with cut-off dis- 



tance rt 



5tq , where = 2 . In reduced MD units 



with LJ parameters a = 1 , e = 1 it is given by 



2\ 4 r -10 -41 

r c ) [r u -r J 



(32) 



Unlike to standard LJ potential function the qLJ potential 
Ea. (l32> has "good" properties at cut-off radius r = r c rj 
2.51 , resulting in better energy conservation. The atom 
mass is assumed to be m = 48 . Number of atoms is fixed 



N = 3375 and all of them are located inside the MD sim- 
ulation cube L = 17 MD units with imposed periodical 
conditions. Before testing the MD system is thermalized 
at T — 1.2 and stored on disk. This thermodynamically 
equilibrium state is a starting point of all testing runs. 

For comparison purposes we choose Verlet algorithm 
in the coordinate form Eqs. (I33I34> : 

Sx h = 6x + a h 2 + 0{h 4 ) (33) 
Xh = x + Sx h (34) 
Oh = f(xh)/m (35) 
v h = (Sx h + 5x Q )/2h + (ah + 2a )V3 + 0(h 4 ) (36) 

To improve the local energy conservation in the Verlet 
scheme the particle velocity Vh is evaluated by using ac- 
curate Ea. J36> which does not affect the trajectory evalua- 
tion (see 1 3)). 

The new integrator is chosen in the form of Eqs.jl5l 
l30l 13 1> with step-by-step algorithm from the Table 2. Un- 
fortunately, the new method requires, in principle, the dou- 
ble run over all i e {1, ..,N} particles to evaluate a^ 
and Ci — Cij . According to our simulations it takes ap- 
proximately 2 (from 1 .8 to 2.2 on different machines) times 
more computer time than it takes for the Verlet scheme. 

Figures 1,2 show the timing data of several simula- 
tion runs for different time steps. Fig.l indicates that root- 
mean-square accuracy of the new method is the superior, 
but the new method takes almost twice as much of CPU 
time. We evaluate the relative efficiency of the method by 
comparing timings required for simulation of the same MD 
time by the new method and Verlet method at given accu- 
racy. The results are demonstrated in Fig. 2 . For given test 
the new integrator becomes more efficient if the desired ac- 
curacy is better than 5 x 10~ 5 . For instance, at the pre- 



scribed accuracy 10 6 the new method takes 4 times less 
computer time than the Verlet algorithm. 

CONCLUSION 

We have demonstrated that the new integrator is time- 
reversible and very accurate. It can be efficiently applied 
for highly precise MD simulations. The new method is ex- 
pected to be useful for simulation of celestial mechanics 
problems too. 

It should be noted that Ea. (l31> was obtained first in 
|4 |, but it was used in a time-irreversible algorithm. 
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